
COMPILE_OPT IDL2, LOGICAL_PREDICATE

COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
COMMON VEL_BLOCK, v, v_std
COMMON HRDS_BLOCK, hrds, hdist


conf = omni_load_conf()
s = omni_read_cat(conf.cat,ncat)
restore,'./local/'+conf.survey+'_velocities.sav',/ver
v_std = findgen(conf.nvbin)*conf.deltav + conf.vstart


dpdfs = omni_read_conffile('./conffiles/dpdf_params.conf')
ancil = omni_read_conffile('./conffiles/ancillary.conf')
d = dindgen(dpdfs.nbins)*dpdfs.binsize + dpdfs.binstart

IF ~exist(pvec) THEN restore,'./local/'+conf.survey+'_pvec.sav',/ver
stat = pvec.stat

t_start = systime(1)

nloop = 5

ndr = fix(floor(4*5*1.5))
ndv = fix(floor(5*5*1.5))


dr = (dindgen(ndr)*1.d + 1.d)*1.d
dv = dindgen(ndv)*0.2d + 2.d

goto,plot_only

blank = lonarr(ndr,ndv)

nhii = blank
ngt1 = blank
ngt1u = blank

fgt1 = double(blank)
fgt1u = double(blank)

nloop = ndv*ndr

FOR ii=0,ndr-1 DO BEGIN
   FOR jj=0,ndv-1 DO BEGIN
      
      ll = jj + ii*ndv
      
      undefine,n_hrds,kdars,pvec
      n_hrds = lonarr(ncat)
      kdars  = bytarr(ncat)
      pvec   = dblarr(ncat,1000)
      
      print,'---------------------------------------'
      print,'Grid points: ',dr[ii],dv[jj]
      
      ;; Set parameter values for this run
      ancil.ppv_dis = dr[ii]
      ancil.ppv_dv  = dv[jj]
      undefine,adist,hdist
      
      
      FOR kk=0,ncat-1 DO BEGIN
         prob = prob_hrds(s[kk],N_HRDS=nhrds,KDAR=kdar)
         n_hrds[kk]  = nhrds
         kdars[kk]   = kdar
         pvec[kk,*]  = prob
      ENDFOR
      
      ;; Number of sources with HRDS association
      ind = where(n_hrds NE 0, ni)
      print,'N w/ HRDS: ',ni
      print,'NHRDS: ',m4_stat(n_hrds[ind])
      nhii[ii,jj] = ni
      
      ;; Number of sources with multiple HRDS association
      ind = where(n_hrds GT 1, ni)
      print,'N mult: ',ni
      ngt1[ii,jj] = ni
      fgt1[ii,jj] = double(ni)/nhii[ii,jj]
      print,'Multiple Fraction: ',fgt1[ii,jj]
      
      IF ni GT 0 THEN BEGIN
         ;; Analyze the KDAR & PVEC
         kd = kdars[ind]        ; KDARs
         p  = pvec[ind,*]       ; PVECs
         st = stat[ind]         ; STATs
         
         ind = where(kd EQ 3 OR kd EQ 7, nmbad)
         ngt1u[ii,jj] = nmbad
         fgt1u[ii,jj] = double(nmbad) / ngt1[ii,jj]
         
         FOR mm=0,ni-1 DO BEGIN
            
            ;; Skip for "pure" KDAR
            IF (kd[mm] EQ 1) || (kd[mm] EQ 2) || (kd[mm] EQ 4) THEN CONTINUE
            
            ;; Analyze the Mutts.
            CASE 1 OF
               ;; N + T OR F + T -- Still removes probability at other KD
               kd[mm] EQ 5 || kd[mm] EQ 6: 
               
               ;; N + F OR N + F + T -- Need to check probability split
               kd[mm] EQ 3 || kd[mm] EQ 7: BEGIN
                  
                  ;; These are defined to be bad, based on looking at
                  ;; the DPDFs.  So, count these?
                  
               END
               ELSE: print,'WTF????'
            ENDCASE
         ENDFOR
      ENDIF
      
      etd = t_start + (systime(1) - t_start) / (double(ll+1) / nloop)
      message,'ETD: '+systime(0,etd)+'  ['+$
              string(ll+1,nloop,format="(I6,' / ',I6)")+']',/inf
      
      undefine,n_hrds,kdars,pvec
   ENDFOR
ENDFOR

save,nhii,ngt1,ngt1u,fgt1,fgt1u,filename='./grs_paper/hrds_mult.sav',/ver

plot_only:

restore,'./grs_paper/hrds_mult.sav',/ver

ddr = arrstep(dr)
ddv = arrstep(dv)

;; Do stuff!

myps,'./grs_paper/hrds_mult.eps',xsize=10,ysize=10
multiplot_xm,[2,2],gap=0.035,mpcharsize=1.0,/doxaxis,/doyaxis

;;==============================
;;  Panel the First

cgPlot,/nodata,[0],xst=4,yst=4,tit='N BGPS w/ Multiple HRDS',charsize=1.0
cgLoadct,41
plotimage,ngt1,range=set_plot_range(ngt1),xst=4,yst=4,min_dpi=100

cgAxis,xaxis=0,xr=minmax(dr)+[-1,1]*ddr/2.,/xst,/save,charsize=1.0,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=minmax(dv)+[-1,1]*ddv/2.,/yst,/save,charsize=1.0,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

cgContour,ngt1,dr,dv,/over,color='black',levels=findgen(6)*100.
help,ngt1,dr,dv
print,m4_stat(dr,dv)

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1


multiplot,/doxaxis,/doyaxis


;;==============================
;;  Panel the Second

cgPlot,/nodata,[0],xst=4,yst=4,tit='N BGPS w/ Multiple HRDS - Unconstrained',$
       charsize=1.0
cgLoadct,41
plotimage,ngt1u,range=set_plot_range(ngt1u),xst=4,yst=4,min_dpi=100

cgAxis,xaxis=0,xr=minmax(dr)+[-1,1]*ddr/2.,/xst,/save,charsize=1.0,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=minmax(dv)+[-1,1]*ddv/2.,/yst,/save,charsize=1.0,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

cgContour,ngt1u,dr,dv,/over,color='black',levels=findgen(7)*20.

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1




multiplot,/doxaxis,/doyaxis

;;==============================
;;  Panel the Third

cgPlot,/nodata,[0],xst=4,yst=4,tit='Fraction of HRDS sample with Multiples',$
       charsize=1.0
cgLoadct,41
plotimage,fgt1,range=set_plot_range(fgt1),xst=4,yst=4,min_dpi=100

cgAxis,xaxis=0,xr=minmax(dr)+[-1,1]*ddr/2.,/xst,/save,charsize=1.0,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=minmax(dv)+[-1,1]*ddv/2.,/yst,/save,charsize=1.0,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

cgContour,fgt1,dr,dv,/over,color='black',levels=findgen(11)*0.05

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1


multiplot,/doxaxis,/doyaxis

;;==============================
;;  Panel the Fourth

cgPlot,/nodata,[0],xst=4,yst=4,tit='Fraction of Panel 1 Unconstrained',$
       charsize=1.0
cgLoadct,41
plotimage,fgt1u,range=set_plot_range(fgt1u),xst=4,yst=4

cgAxis,xaxis=0,xr=minmax(dr)+[-1,1]*ddr/2.,/xst,/save,charsize=1.0,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=minmax(dv)+[-1,1]*ddv/2.,/yst,/save,charsize=1.0,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

cgContour,fgt1u,dr,dv,/over,color='black',levels=findgen(11)*0.05

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1




myps,/done,/mp

END
